rm(list=ls())
setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")
library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)

icb<-read.csv("icb2v13.csv")
icb<-subset(icb, icb$systrgyr>1945)
icb_nuc_stakes<-icb[,c("crisno", "systrgyr", "cracid",
                       "actor", "nuclear", "gravty", "outcom")]

colnames(icb_nuc_stakes)<-c("crisno", "year", "statea",
                            "namea", "nuclear_a", "gravty_a", "outcom_a")
icb_nuc_stakes$namea<-ifelse(icb_nuc_stakes$statea==484, "CON", icb_nuc_stakes$namea)
icb_nuc_stakes$namea<-ifelse(icb_nuc_stakes$statea==490, "DRC", icb_nuc_stakes$namea)


icbdy<-read.csv("icbdy_v13.csv")
icbdy<-subset(icbdy, icbdy$year>1945 & icbdy$ongoing==0)
icbdy<-icbdy[,c("crisno", "crisname", 
                "statea", "namea",
                "stateb", "nameb")]
#install.packages("countrycode")
library(countrycode)
icbdy$namea<-countrycode(icbdy$statea, "cown", "cowc")
icbdy$namea<-ifelse(icbdy$statea==260, "GFR", icbdy$namea)
icbdy$nameb<-countrycode(icbdy$stateb, "cown", "cowc")

## Changing the Data to reflect that the DRC, not Congo-Brazzaville, was involved 
# in the Congo II crisis
icbdy$stateb<-ifelse(icbdy$crisno==211 & icbdy$stateb==484, 490, icbdy$stateb)
icbdy$nameb<-ifelse(icbdy$crisno==211 & icbdy$stateb==490, "DRC", icbdy$nameb)

test<-merge(icbdy, icb_nuc_stakes, by=c("crisno", "statea", "namea"),
            all.x = T)

colnames(icb_nuc_stakes)<-c("crisno", "year", "stateb",
                            "nameb", "nuclear_b", "gravty_b", "outcom_b")
test2<-merge(test, icb_nuc_stakes, by=c("crisno", "stateb", "nameb"),
            all.x = T)

nuc_cris<-subset(test2, test2$nuclear_a>2 | test2$nuclear_b>2)

nuc_cris$outcom_a<-ifelse(is.na(nuc_cris$outcom_a)==T, 0, nuc_cris$outcom_a)
nuc_cris$outcom_b<-ifelse(is.na(nuc_cris$outcom_b)==T, 0, nuc_cris$outcom_b)
nuc_cris$gravty_a<-ifelse(is.na(nuc_cris$gravty_a)==T, 0, nuc_cris$gravty_a)
nuc_cris$gravty_b<-ifelse(is.na(nuc_cris$gravty_b)==T, 0, nuc_cris$gravty_b)
nuc_cris$nuclear_a<-ifelse(is.na(nuc_cris$nuclear_a)==T & nuc_cris$statea==2, 4, nuc_cris$nuclear_a)
nuc_cris$nuclear_a<-ifelse(is.na(nuc_cris$nuclear_a)==T, 1, nuc_cris$nuclear_a)
nuc_cris$nuclear_b<-ifelse(is.na(nuc_cris$nuclear_b)==T & 
                           nuc_cris$stateb==365 & nuc_cris$year.x<1950, 2, nuc_cris$nuclear_b)
nuc_cris$nuclear_b<-ifelse(is.na(nuc_cris$nuclear_b)==T & 
                             nuc_cris$stateb==365 & nuc_cris$year.x>1958, 4, nuc_cris$nuclear_b)
nuc_cris$nuclear_b<-ifelse(is.na(nuc_cris$nuclear_b)==T, 1, nuc_cris$nuclear_b)


nuc_cris$gravty_b<-ifelse(is.na(nuc_cris$gravty_b)==T, 0, nuc_cris$gravty_b)

library(dplyr)
nuc_cris<- nuc_cris %>% group_by(statea, stateb) %>% mutate(dyid=paste(max(statea, stateb), min(statea, stateb),sep='_'))
nuc_cris<- nuc_cris %>% group_by(dyid, crisno) %>% mutate(crdyid=paste(dyid, crisno, sep = '_'))

# Deleting the second Yom Kippur Observation, where Israel wins, since that is what Kroenig did
nuc_cris$second_yom<-ifelse(nuc_cris$crisno==255 & nuc_cris$stateb==666 & nuc_cris$outcom_b==1, 1, 0)
nuc_cris<-subset(nuc_cris, nuc_cris$second_yom==0)

# Deleting the second War of Attrition Observation, where Israel wins, since that is what Kroenig did
nuc_cris$second_woa<-ifelse(nuc_cris$crisno==232 & nuc_cris$stateb==666 & nuc_cris$outcom_b==1, 1, 0)
nuc_cris<-subset(nuc_cris, nuc_cris$second_woa==0)

## Checking sum of nuclear weapons variable
nuc_cris<- nuc_cris %>% group_by(crdyid) %>% mutate(nucdysum=sum(nuclear_a, nuclear_b))

## Eliminating crises between two nuclear states, but ensuring only crises where at least
# one side has nuclear weapons. 
## ICB ranks nuclear capability of states on a four point scale: 3 indicates nuclear
# capability but no second strike, 4 indicates nuclear capability with second strike
## Thus, the dyadic sum of nuclear must be at most six. All of the dyads with six must be
# a combo of a 4 and a 2. 
## The dyadic sum of nuclear must be greater than 3, but all of the fours must be a combo
## of 3 and 1, and not 2 and 2. 
nuc_cris<-subset(nuc_cris, nuc_cris$nucdysum<=6)
nuc_cris<-subset(nuc_cris, nuc_cris$nucdysum>3)
nuc_cris$twothrees<-ifelse(nuc_cris$nucdysum==6 & nuc_cris$nuclear_a==3 & nuc_cris$nuclear_b==3, 1, 0)
nuc_cris<-subset(nuc_cris, nuc_cris$twothrees==0)

## Stakes
nuc_cris$highstakes_a<-ifelse(nuc_cris$gravty_a==2 | nuc_cris$gravty_a==3 | 
                                nuc_cris$gravty_a>4, 1, 0)
nuc_cris$highstakes_b<-ifelse(nuc_cris$gravty_b==2 | nuc_cris$gravty_b==3 | 
                                nuc_cris$gravty_b>4, 1, 0)
nuc_cris$inferior_hs_1<-ifelse(nuc_cris$nuclear_a<3 & nuc_cris$highstakes_a==1 | 
                               nuc_cris$nuclear_b<3 & nuc_cris$highstakes_b==1, 1, 0)
nuc_cris$nuclear_hs_1<-ifelse(nuc_cris$nuclear_a>2 & nuc_cris$highstakes_a==1 | 
                               nuc_cris$nuclear_b>2 & nuc_cris$highstakes_b==1, 1, 0)
nuc_cris$bothhighstakes<-ifelse(nuc_cris$inferior_hs_1==1 & nuc_cris$nuclear_hs_1==1, 1, 0)
nuc_cris$inferior_hs<-ifelse(nuc_cris$inferior_hs_1==1 & nuc_cris$nuclear_hs_1==0, 1, 0)
nuc_cris$nuclear_hs<-ifelse(nuc_cris$inferior_hs_1==0 & nuc_cris$nuclear_hs_1==1, 1, 0)
nuc_cris$nohighstakes<-ifelse(nuc_cris$inferior_hs_1==0 & nuc_cris$nuclear_hs_1==0, 1, 0)

## Table 2 from Appendix I 
nuc_cris<- nuc_cris %>% group_by(crdyid) %>% mutate(outcom_sum=sum(outcom_a, outcom_b))

nuclear_outcomes<-table(nuc_cris$outcom_a[nuc_cris$nuclear_a>2])+table(nuc_cris$outcom_b[nuc_cris$nuclear_b>2])
non_nuclear_outcomes<-table(nuc_cris$outcom_a[nuc_cris$nuclear_a<3])+table(nuc_cris$outcom_b[nuc_cris$nuclear_b<3])

outcome_mat<-matrix(NA, 3, 5)
rownames(outcome_mat)<-c("nuclear outcomes", "non nuclear outcomes", "diff in means")
colnames(outcome_mat)<-c("No Outcome Coded", "Victory", "Compromise", "Stalemate", "Defeat")
outcome_mat[1,]<-c(0,nuclear_outcomes)/155
outcome_mat[2,]<-non_nuclear_outcomes/155
outcome_mat[3,]<-outcome_mat[1,]-outcome_mat[2,]

xtable(round(outcome_mat, digits=2))



## Table 3 from Appendix I
nuc_outcomes_bhs<-table(nuc_cris$outcom_a[nuc_cris$bothhighstakes==1 & nuc_cris$nuclear_a>2])+
  c(table(nuc_cris$outcom_b[nuc_cris$bothhighstakes==1 & nuc_cris$nuclear_b>2]),0)

nuc_outcomes_nhs<-table(nuc_cris$outcom_a[nuc_cris$nuclear_hs==1 & nuc_cris$nuclear_a>2])+
  table(nuc_cris$outcom_b[nuc_cris$nuclear_hs==1 & nuc_cris$nuclear_b>2])

nuc_outcomes_ihs<-table(nuc_cris$outcom_a[nuc_cris$inferior_hs==1 & nuc_cris$nuclear_a>2])+
  c(table(nuc_cris$outcom_b[nuc_cris$inferior_hs==1 & nuc_cris$nuclear_b>2]),0)

nuc_outcomes_nohs<-table(nuc_cris$outcom_a[nuc_cris$nohighstakes==1 & nuc_cris$nuclear_a>2])+
  table(nuc_cris$outcom_b[nuc_cris$nohighstakes==1 & nuc_cris$nuclear_b>2])

nuc_outcomes_stakes<-matrix(NA, 4, 5)
rownames(nuc_outcomes_stakes)<-c("Both High Stakes", "Only Nuclear State High Stakes",
                                  "Only Non-Nuclear State High Stakes", "Neither High Stakes")
colnames(nuc_outcomes_stakes)<-c("No Outcome Coded", "Victory", "Compromise", "Stalemate", "Defeat")
nuc_outcomes_stakes[1,]<-c(0,nuc_outcomes_bhs)/26
nuc_outcomes_stakes[2,]<-c(0, nuc_outcomes_nhs,0)/15
nuc_outcomes_stakes[3,]<-c(0, nuc_outcomes_ihs)/70
nuc_outcomes_stakes[4,]<-c(0, nuc_outcomes_nohs)/44
nuc_outcomes_stakes

xtable(round(nuc_outcomes_stakes, digits=2))

## Table 4 fro mAppendix I

non_nuc_outcomes_bhs<-c(0,table(nuc_cris$outcom_a[nuc_cris$bothhighstakes==1 & nuc_cris$nuclear_a<3]))+
 c(0, table(nuc_cris$outcom_b[nuc_cris$bothhighstakes==1 & nuc_cris$nuclear_b<3]))

table(nuc_cris$outcom_a[nuc_cris$nuclear_hs==1 & nuc_cris$nuclear_a<3])
table(nuc_cris$outcom_b[nuc_cris$nuclear_hs==1 & nuc_cris$nuclear_b<3])
tab_replace<-matrix(NA,1,5)
tab_replace<-c(7,2,1,4,1)
non_nuc_outcomes_nhs<-tab_replace

non_nuc_outcomes_ihs<-c(0,table(nuc_cris$outcom_a[nuc_cris$inferior_hs==1 & nuc_cris$nuclear_a<3]))+
  c(0,table(nuc_cris$outcom_b[nuc_cris$inferior_hs==1 & nuc_cris$nuclear_b<3]))

non_nuc_outcomes_nohs<-table(nuc_cris$outcom_a[nuc_cris$nohighstakes==1 & nuc_cris$nuclear_a<3])+
  table(nuc_cris$outcom_b[nuc_cris$nohighstakes==1 & nuc_cris$nuclear_b<3])

non_nuc_outcomes_stakes<-matrix(NA, 4, 5)
rownames(non_nuc_outcomes_stakes)<-c("Both High Stakes", "Only Nuclear State High Stakes",
                                 "Only Non-Nuclear State High Stakes", "Neither High Stakes")
colnames(non_nuc_outcomes_stakes)<-c("No Outcome Coded", "Victory", "Compromise", "Stalemate", "Defeat")
non_nuc_outcomes_stakes[1,]<-non_nuc_outcomes_bhs/26
non_nuc_outcomes_stakes[2,]<-non_nuc_outcomes_nhs/15
non_nuc_outcomes_stakes[3,]<-non_nuc_outcomes_ihs/70
non_nuc_outcomes_stakes[4,]<-non_nuc_outcomes_nohs/44
non_nuc_outcomes_stakes

xtable(round(non_nuc_outcomes_stakes, digits=2))


## Table 5 from Appendix I
stakes_outcome_table<-matrix(NA, 4, 5)
rownames(stakes_outcome_table)<-c("Both High Stakes", "Only Nuclear State High Stakes",
                                  "Only Non-Nuclear State High Stakes", "Neither High Stakes")
colnames(stakes_outcome_table)<-c("No Outcome Coded", "Victory", "Compromise", "Stalemate", "Defeat")
stakes_outcome_table[1,]<-(c(0,nuc_outcomes_bhs)/26)-(non_nuc_outcomes_bhs/26)
stakes_outcome_table[2,]<-((c(0,nuc_outcomes_nhs,0))/15)-(non_nuc_outcomes_nhs/15)
stakes_outcome_table[3,]<-(c(0,nuc_outcomes_ihs)/70)-(non_nuc_outcomes_ihs/70)
stakes_outcome_table[4,]<-(c(0,nuc_outcomes_nohs)/44)-(non_nuc_outcomes_nohs/44)
stakes_outcome_table


xtable(round(stakes_outcome_table, digits=2))
